Saturation of Elliptic Flow and the Transport Opacity 
of the Gluon Plasma at RHIC 



Denes Molnar ' 2 and Miklos Gyulassy 1 ' 2,3 
1 Physics Department, Columbia University, 538 W. 120th Street, New York, NY 10027 
2 RMKI Research Institute for Particle and Nuclear Physics, PO Box 49, H-1525 Budapest, Hungary 
3 Collegium Budapest, Szentharomsag u. 2., H-1014 Budapest, Hungary 
(Final published version + NPA703 Erratum) 



o 
o 

(N 
o 

Q 

m 

m 
> 
m 
o 
o 
^1- 
o 



o 

=5 



X 



Differential elliptic flow and particle spectra are calculated 
taking into account the finite transport opacity of the gluon 
plasma produced in Au+Au at E cm ~ 130 A GeV at RHIC. 
Covariant numerical solutions of the ultrarelativistic Boltz- 
mann equation are obtained using the MPC parton cascade 
technique. For typical pQCD (~ 3 mb) elastic cross sections, 
extreme initial gluon densities, dN g /dn ~ 15000, are required 
to reproduce the elliptic flow saturation pattern reported by 
STAR. However, we show that the solutions depend mainly 
on the transport opacity, X ~ J dzatpg, and thus the data can 
also be reproduced with dN g /dr] ~ 1000, but with extreme 
elastic parton cross sections, ~ 45 mb. We demonstrate that 
the spectra and elliptic flow are dominated by numerical ar- 
tifacts unless parton subdivisions ~ 100 — 1000 are applied to 
retain Lorentz covariance for RHIC initial conditions. 

PACS 12.38. Mh; 24.85. +p; 25. 75. - q 



I. INTRODUCTION 

Differential elliptic flow, V2(p±) = (cos(2<f>)) p± , the sec- 
ond Fourier moment of the azimuthal momentum distri- 
bution for fixed p± , is one of the important experimental 
probes of collective dynamics in A + A reactions [1-13]. 
The discovery [14] of a factor of two enhancement of el- 
liptic flow in noncentral nuclear collisions at RHIC rela- 
tive to SPS [15] has generated even more interest in this 
"barometric" measure of collective transverse flow. In ad- 
dition, preliminary STAR data reported in [16] suggest 
a remarkable saturation property of this flow at high p± 
with V2(p± > 2 GeV) ~ 0.15. This corresponds to a fac- 
tor of two azimuthal angle asymmetry of high-pj_ particle 
production relative to the reaction plane. This collec- 
tive effect depends strongly on the dynamics in a heavy 
ion collision and therefore provides important constraints 
about the density and effective energy loss of partons. In 
particular, we show that it constrains the transport opac- 
ity of the produced gluon plasma. 

Predictions of collective elliptic flow in noncentral nu- 
clear collisions were first based on ideal (nondissipative) 
hydrodynamics [1,2,6]. Unlike at lower energies, ideal hy- 
drodynamics seems to reproduce the (low px < 2 GeV) 
data [14] at RHIC remarkably well. However, it fails to 
saturate at high p± > 2 GeV as indicated by the prelim- 
inary data [16]. The hydrodynamic results were found 
in [6] to be surprisingly insensitive to the choice of ini- 



tial conditions, equation of state and freezeout criteria, 
once the observed dN c h/drj was reproduced, leaving no 
adjustable hydrodynamic model parameters with which 
the saturation property could be reproduced. 

The lack of saturation in ideal hydrodynamics is due 
to the assumption of zero mean free path and that local 
equilibrium can be maintained until a chosen freezeout 
3D hypersurface is reached. This idealization is certainly 
invalid outside some finite domain of phase space in heavy 
ion collisions [17]. Finite transition rates are expected 
to produce nonequilibrium deviations from the predicted 
hydrodynamic flow pattern. Covariant Boltzmann trans- 
port theory provides a convenient framework to estimate 
dissipative effects. The assumption of local equilibrium 
is replaced by the assumption of a finite local mean free 
path X(s,x) = l/a(s)n(x). The theory then naturally 
interpolates between free streaming (A = oo) and ideal 
hydrodynamics (A = 0). 

The previous calculations of collective flow from co- 
variant transport theory [9,18,19] lead to too small col- 
lective effects. This was due to the use of small per- 
turbative QCD cross sections and dilute parton initial 
densities based on HIJING [20] . Recently, denser parton 
initial conditions were suggested based on gluon satura- 
tion models [21]. Initial gluon densities up to five times 
higher than from HIJING were predicted. The question 
studied in this paper is whether such initial conditions 
may be dense enough to generate the observed collective 
flow even with pQCD elastic cross sections. In this pa- 
per, we explore the dependence of differential elliptic flow 
on initial conditions, or equivalently 1 , on the magnitude 
of partonic cross section. 

We note that most hadronic cascade models supple- 
mented with string dynamics [4,11] underpredict ellip- 
tic flow because the spatial asymmetry is too small after 
hadronization to generate the observed momentum space 
asymmetry. The reason is that hadronization through 
longitudinal string excitations reduces the strength of 
partonic level elliptic flow. To produce sufficient ellip- 
tic flow in string models, other mechanisms have to be 
included such as color exchange. While the impact pa- 
rameter dependence of elliptic flow at RHIC can be repro- 



x The equivalence is due to the scaling property explained in 
Section II B. 
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duced via color exchange [10], it is not known whether the 
saturation property studied here can also be explained. 

The saturation and eventual decrease of t>2 at high p± 
has been predicted as a consequence of finite inelastic 
parton energy loss [12,13]. Those predictions however 
assumed the validity of an Eikonal approach at moder- 
ate pt ~ 10 GeV. In addition, the pQCD computable jet 
quenched part had to be joined phenomenologically onto 
a parametrized "soft" nonperturbative component below 
Pt < 2 GeV/c. Covariant transport theory overcomes 
the need to treat soft and hard dynamics on different 
footings. It is the only practical self-consistent theoreti- 
cal tool at present to address simultaneously both the soft 
collective component and the far from equilibrium high- 
p± component. While current parton cascade techniques 
lack at present a practical means to implement covariant 
inelastic energy loss, it is of considerable theoretical in- 
terest to solve covariant Boltzmann theory even in the 
elastic limit since so few solutions are known. We solve 
that theory numerically here to get insight into the dy- 
namical interplay between the soft and hard components 
over a wide dynamical range of parameters. 

For large enough elastic opacities, the observed collec- 
tive flow strength can certainly be reproduced [19]. The 
outstanding question which we focus on is whether the 
detailed pattern of deviations from ideal hydrodynamic 
flow and its saturation at high p± can also be understood 
quantitatively in this particular dynamical framework. 

Forerunners of this study [9,19] computed elliptic flow 
for partonic systems starting from initial conditions ex- 
pected at RHIC. In this paper we extend those studies 
in three aspects. We compute the p^-diffcrcntial ellip- 
tic flow v 2 {p\_)- The consequences of hadronization are 
investigated, which is necessary to compare to the obser- 
vations. Finally, we use realistic diffuse nuclear geometry 
for the initial conditions. 

We compute the partonic evolution with MPC [22], a 
newly formulated, covariant, parton kinetic theory tech- 
nique. MPC is an extension of the covariant parton cas- 
cade algorithm, ZPC [23] . Both MPC and ZPC have been 
extensively tested [24,25] and compared to analytic trans- 
port solutions and covariant Euler and Navier-Stokes dy- 
namics in 1+1D geometry. A critical feature of both 
these algorithms is the implementation of the parton sub- 
division technique proposed by Pang [25,26]. 

Extensions of MPC to include inelastic 2^-3 partonic 
processes [18] are under development. In this paper, we 
apply MPC in the pure elastic parton interactions mode 
as in ZPC [27]. 



II. COVARIANT TRANSPORT THEORY 



transport theory in which the on-shell phase space den- 
sity f(x, p), evolves with an elastic 2 — > 2 rate as 

Pld^fl = fffifofa ~ hfo) W / 12^34<5 4 (pi +p 2 - P3 - Pi) 
234 

+ S(x,pi). (1) 



Here W is the square of the scattering matrix element, 
the integrals are shorthands for J = J , where 

i 

g is the number of internal degrees of freedom, while 
fj = f( x ,Pj)- The initial conditions are specified by the 
source function S(x,p), which we discuss at the end of 
this Subsection. 

For our applications below, we interpret f(x,p) as 
describing an ultrarelativistic massless gluon gas with 
g = 16 (8 colors, 2 helicities). We neglect quark degrees 
of freedom because at RHIC gluons are more abundant. 

In principle, the transport equation (1) could be ex- 
tended for bosons with the substitution /i/ 2 — > /i/2(l + 
/3)(1 + fi) an d a similar one for 73/4. In practice, no co- 
variant algorithm yet exists to handle such nonlincaritics. 
We therefore limit our study to quadratic dependence of 
the collision integral on /. 

The elastic gluon scattering matrix elements in dense 
parton systems are modeled by a Debye-screened form as 
in Ref. [9]: 



dt 



(2) 



where [i is the screening mass, o~o(s) = 97ra 2 (s)/2/z 2 is 
the total cross section. For simplicity, we will assume o~o 
to be energy independent, neglecting its weak logarithmic 
dependence on s in the energy range relevant at RHIC. 

For small values of fi, forward-peaked scattering is fa- 
vored, while as fi increases the angular distribution be- 
comes more and more isotropic. For a fixed total cross 
section 2 , the relevant transport cross section is 



(3) 



0t(s) = J da et sin 2 6» cm = J dt - 

= a 4z(l + z)[(2z + l)ln(l + l/z)-2] , 

where z = H 2 /s. This is a monotonic function of n and 
maximal in the isotropic (/1 — > 00) case. In the small 
angle dominated limit, with z« 1, <Jt/&o ~ 4z(lnl/z — 
2). 

It is important to emphasize that while the cross sec- 
tion suggests a geometrical picture of action over fi- 
nite distances, we use Eq. (2) only as a convenient 
parametrization to describe the effective local transition 



A. Transport equation 



We consider here, as in Refs. [26,23,22,17], the sim- 
plest but nonlinear form of Lorentz-covariant Boltzmann 



2 To keep the total cross section constant as a function of fi, 
one of course has to adjust the coupling a s accordingly. 
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probability, W. In the present study this is simply mod- 
eled as W(t) — s da/dt. The particle subdivision tech- 
nique (see next Subsection) needed to recover covariance 
removes all notion of nonlocality in this approach, just 
like in hydrodynamics. Thus, the cross sections, e.g., 100 
mb, used in the present study to simulate a high collision 
rate do not imply acausal action at a distance. 

For Au+Au collisions, the initial condition was taken 
to be a longitudinally boost invariant Bjorken tube in 
local thermal equilibrium at temperature To at proper 
time t = 0.1 fm/c with uniform pseudorapidity r\ = 
l/21og((f + z)/(t — z)) distribution between \rf\ < 5. 
The transverse density distribution was assumed to be 
proportional to the binary collision distribution for two 
Woods-Saxon distributions. For collisions at impact pa- 
rameter b the transverse binary collision profile is 



dN(b) ( 

= a jet T A ^ ± 



(4) 



where T^b) = J dzpA(Vz 2 + b 2 ), in terms of the dif- 
fuse nuclear density pa{^)- The pQCD jet cross section 
normalization, <Jj e t, and the temperature Tq were deter- 
mined by fitting the gluon minijet transverse momentum 
spectrum predicted by HIJING [20] for central Au+Au 
collisions at y/s = 130 A GeV (without shadowing and jet 
quenching). This gives dN(Q)/dr} = 210 and T = 700 
McV. 

Evolutions from different initial densities (but the same 
density profile) can be obtained by varying the cross sec- 
tion only and using the scaling property explained in Sec- 
tion II C. 



B. Parton Subdivision 

We utilize the parton cascade method to solve the 
Boltzmann transport equation (1). A critical drawback 
of all cascade algorithms is that they inevitably lead to 
numerical artifacts that violate Lorentz covariance. This 
occurs because particle interactions are simulated to oc- 
cur whenever the distance of closest approach (in the 
relative cm.) is d < ^Jos/it. 

Acausal (superluminal) propagation due to action at a 
distance leads to severe numerical artifacts. In particular, 
the transverse energy evolution dE T (T)/dy and the final 
asymptotic transverse energy per unit rapidity are frame 
dependent [17]. 

To recover the local character of equation (1) and hence 
Lorentz covariance, it is essential to use the parton sub- 
division technique [26,23]. This technique is based on the 
covariance of Eq. (1) under the transformation 

/-/' = */, W^W' = W/£ (a - a' = a/t) , (5) 

where I is the number of particle subdivisions. The mag- 
nitude of numerical artifacts is governed by the dilutc- 
ness of the system y/a /Xmfp, which scales with \j\fl 



[25] . Lorentz violation therefore formally vanishes in the 
I — > 00 limit. The convergence to the accurate covariant 
solution with I is slower if the density or cross section 
increases. 

Figures 1 and 2 illustrate the effect of Lorentz violation 
on the observables studied in this paper. For insufficient 
particle subdivision, elliptic flow and the p± spectra are 
dominated by numerical artifacts. In particular, ellip- 
tic flow is significantly underpredicted, while the high- 
p± spectrum exhibits an unphysical "reheating" during 
the expansion. For the initial conditions for these plots, 
these numerical artifacts disappear only when the parti- 
cle subdivision reaches ~ 200. This reinforces the results 
by Ref. [17], where very high ~ 100 — 1000 subdivisions 
were needed to obtain accurate numerical solutions of the 
transport equation (1) for initial conditions expected at 
RHIC. 



C. Scaling of the transport solutions 

Subdivision covariance (5) actually implies that the 
transport equation has a broad dynamical range, and 
the solution for any given initial condition and transport 
property immediately provides the solution to a broad 
band of suitably scaled initial conditions and transport 
properties. This is because solutions for problems with 
I times larger the initial density dN/di]d 2 x±, but with 
one £-th the reaction rate can be mapped to the original 
(£ = 1) case for any I. We must use subdivision to elimi- 
nate numerical artifacts. However, once that is achieved, 
we have actually found the solution to a whole class of 
suitably rescaled problems. 

The dynamical range of the transport equation (1) is 
further increased by its covariance under coordinate and 
momentum rescaling [17], leading to covariance of the 
transport theory under 



m — > ml = m/£p, 



(6) 



where £ x and l v are the coordinate and momentum scal- 
ing parameters, respectively. This means [17] that we 
can scale one solution to others provided that both p/To 
and aodN/dij remain the same (we cannot exploit co- 
ordinate scaling because the nuclear geometry is fixed). 
For example, three times the density with one-third the 
cross section leaves both parameters the same, hence the 
results can be obtained via scaling without further com- 
putation. 

In general the numerical (cascade) solution of Eq. (1) 
tends in the I — > 00 limit toward a covariant physical so- 
lution that depends on two scales, h/Tq and aodN/dr). In 
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an Eikonal picture of high-pj_ production, the distribu- 
tions are expected to depend on the opacity or the mean 
number of collisions in the medium 



(n) 



L 

dN a Rg 
dy 2nR 2 G ° g t 



(7) 



where r is the formation proper time and Rq is the 
effective Gaussian transverse coordinate rms radius. 

Our numerical results in Table I show that for a given 
centrality and initial gluon density, the average number of 
collisions per parton is within 10% accuracy proportional 
to do an d does not depend on /i/Tq: 

. dN(0) n , x a Q , n . dN(0) n! NX 



dr] T 



Therefore, one would naively expect a^dN/dn to be the 
relevant scale. 

However, from the point of view of dissipative dynam- 
ics via Navier-Stokes and Fokker-Planck equation, the 
more relevant dynamical parameter is the effective trans- 
port opacity 



Of / \ 



( J dzp(x + zn,T=^j) . (9) 



The ensemble average over initial coordinates and direc- 
tions is implied above. In addition, o~ t here stands for 
o~l = (Tt((s)), where (s) — I8T5 2 is the initial thermal 
average of s. This is an approximation to the ensemble 
average (o- t (s)). 

In general, the transport opacity is a dynamical quan- 
tity that we do not know until we solved the transport 
equation for the set of parameters b, o~odN(0)/dn, and 
o~t/o~o ( or equivalently, /J, /To). However, Eqs. (8) and (9) 
imply that for the range of parameters considered in this 
study it is approximately proportional to the product of 
the two scales aodN(0)/dn and o- t /o- : 



dN(0) a t 

X(b,o- —j , — ) 

an a 



; dN'(y) ^ a o^—)) • (io) 

"0 dr] 



The nontrivial, impact parameter dependent part 
(n(b, o~odN(Q) / drj) is the average number of collisions per 
parton as a function of b for a fixed value of crodA(O) /dr) 
(and an arbitrary /jl/Tq), which is tabulated in Table 
I. For example, one could use the values of set E), in 
which case the corresponding proportionality constant is 
a' dN'(0)/dn = 3 mb x 210 = 630 mb. 

Of course, there is no a priori guarantee that the so- 
lutions of Eq. (1) depend only on this transport opacity 
parameter. However, as demonstrated in Fig. 3, this 
turns out to hold within ~ 10 — 20% accuracy for ellip- 
tic flow and the transverse momentum spectra out to 6 
GeV/c for the parameters and initial conditions appro- 
priate at RHIC energies. In particular, simulations with 



very different impact parameters give the same ratio be- 
tween the initial and final spectrum, provided \ is the 
same. Hence, the relevant parameter that governs the 
evolution of the p± spectra is \ alone. This is not so 
for elliptic flow, which is also driven by the initial spatial 
anisotropy and thus depends not only on \ but also on 
the impact parameter. 



III. NUMERICAL RESULTS FOR THE 
PARTONIC EVOLUTION 

In this Section we present elliptic flow results and p± 
spectra for the partonic evolution. 

Under the scaling (6) the differential elliptic flow 
1*2 (pj_) and the p± spectrum transform as 



V2(p±) —> v' 2 (p±) = v 2 



P-L 



dN 
d 2 p± 



(P±) 



dN ' 
d 2 p± 



(P±) 



, dN 
d 2 p± 



P-L 



^ • (11) 



Hence elliptic flow depends on a t and dN/dn only 
through the product a t dN/dr). On the other hand, the 
p± spectrum depends on a t and dN/dn separately. 

As emphasized in the previous Section, we will label 
the results by the effective elastic transport opacity \ 
from Eq. (9) and the impact parameter b. Possible initial 
gluon densities and transport cross sections correspond- 
ing to a given transport opacity \ and impact parameter 
b can be extracted using Eq. (10), while possible cutoffs 
and total cross sections corresponding to a given trans- 
port cross section can be obtained from Eq. (3). These 
mappings are not unique. A given x and b correspond 
to a whole class of possible initial densities, total cross 
sections and cutoffs. 

Tabic I shows the set of simulation parameters for each 
simulation, together with \ determined directly from the 
average number of cascade collisions per particle. In Ta- 
ble I we also introduced letter codes A) through F) as 
a quick reference to particular subsets of simulation pa- 
rameters. We will include this letter code on most labels 
together with x, for convenience. 

The evolution was performed numerically with 40 and 
100 mb isotropic cross sections, and with 3, 40 and 100 
mb gluonic cross sections with /i/Tq = 1. We used parti- 
cle subdivision t = 100 for impact parameters 0, 2, and 
4 fm, while I = 220, 450, 1100, and 5000, for 6 = 6, 
8, 10, and 12 fm. Our study of subdivision convergence 
shown in Figs. 1 and 2 indicates that with b = 8 fm 
and x = 9. 74^ a particle subdivision of I ~ 200 - 250 
is required, which means that for b = 0, 2, and 4 fm, 
I = 100 is not sufficient if x > 8 — 10. Unfortunately, our 
computational resources were insufficient to allow higher 
subdivision runs. While this affects elliptic flow results 
little because most elliptic flow contributions come from 
b > 4 fm, the particle spectra are affected significantly. 
Therefore, we only present spectra for b > 6 fm. 
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A. Elliptic flow results 

Figures 4 and 5 show the final asymptotic gluon elliptic 
flow as a function of transverse momentum for different 
impact parameters with a t dN(0)/di] w 2580°' mb and 
g440 A ) mb, respectively (the corresponding transport 
opacities at b = are Xb=o = 7.90 c) and 19.4 A '). With 
increasing elliptic flow increases until p± ~ 1.5 — 2 
GcV, where it saturates, reproducing the pattern ob- 
served at RHIC [14,16]. With increasing impact parame- 
ter, elliptic flow first monotonically increases, then mono- 
tonically decreases, showing a maximum at b w 8 fm. 
These features were universal for all the cross sections 
we studied, except for a small increase in the location of 
the maximum with increasing transport cross section, as 
can be seen in Fig. 6, from b = 7 fm at a t — 0.91 B - ) 
mb to 9 fm at 66 s ) mb. Also, as expected, elliptic flow 
is a monotonically increasing function of the transport 
opacity, if the impact parameter is kept fixed. 

Figure 6 shows the p_L-integrated gluon elliptic flow 
as a function of centrality. The cascade reproduces the 
trend seen in the STAR data [14] down to very small 
centralities ~ 0.1 — 0.2, where the ideal hydrodynamical 
assumption of zero mean free path certainly breaks down. 
To quantitatively reproduce the data, transport opacities 
Xb=o ^8 — 14 are needed. With the pQCD elastic gg 
cross section <7o(/U = T) w 3 mb, this corresponds to an 
initial gluon density dN g (0)/dr) ~ 3000 - 5000. 

Figure 7 shows the impact-parameter-averaged gluon 
elliptic flow as a function of transverse momentum for 
different transport opacities. The impact-parameter- 
averaged flow was computed via the formula 



2tt 



^max JO 



dbb 



/d0 008(20)3^.(6) 



dN 
dr]d 2 p± 



dN_ 

Pi 

(b) 



(12) 



with b max = 12 fm. As we show in the Appendix, for 
our transport theory solutions, Eq. (12) gives compara- 
ble results to the minimum-bias differential elliptic flow 
defined by STAR as 



„STAR 



(P±) 



lo^bdbjd^^ib) 



(13) 



which weights flow in more central events preferentially. 
This is not the case for ideal hydrodynamic solutions [6] , 
for which the STAR definition [14] results in much smaller 
flow than if Eq. (12) is used. 

We use the definition (12) because it can be computed 
numerically more reliably in our approach. As discussed 
in the beginning of Section III, it is difficult to reduce 
cascade numerical artifacts to an acceptable level when 
the transport opacity \ isi large. For all other parameters 
kept fixed, x increases with decreasing impact parameter, 
hence the STAR definition that weights V2 at small b 
preferentially is more prone to such numerical artifacts. 



Varying the magnitude of energy loss we searched for 
the drop in v 2 {p±) at high pj_ predicted by calculations 
based on inelastic parton energy loss [12,13]. Although 
those studies consider only effects due to radiative energy 
loss, one expects a similar behavior in case of purely elas- 
tic energy loss. 

Figure 8 shows the dependence of vi (p± ) on the trans- 
port opacity for a fixed impact parameter. We varied the 
opacity by changing the screening mass [i. As expected, 
elliptic flow decreases with decreasing x- However, there 
is no sign of a drop at high p^: within statistical errors, 
the results are consistent with a constant flow from 2 to 
6 GcV transverse momentum. 



B. Particle spectra 

Figures 9-11 show the final gluon p± spectra from MPC 
as a function of transport opacity and impact parameter. 
To show more clearly the degree of quenching due to 
multiple elastic collisions, we plot in Figs. 9 and 10 the 
ratio of the final spectra to the initial for b = 6, 8, 10, and 
12 fm. With HIJING initial densities and pQCD elastic 
cross sections, \ ~ 0.02 — 0.2 is too small to produce more 
than ~ 10% quenching. As we increase the transport 
opacity, quenching of the p± > 2 GeV/c range increases 
and by x ^ 14 — 16 it reaches a factor of ten suppression 
at p ± > 6 GeV. 

While the quenching depends on x only, the absolute 
yield is proportional to the initial dN/dn. Hence from 
the absolutely normalized measured spectrum at a given 
centrality one could extract both cr t and dN(0)/dr). As 
seen in Fig. 11, the quenching at high p± is comple- 
mented by an enhancement at low p±. For clarity, we 
normalized all curves at p± = 2 GeV. While the STAR 
data [30] are too preliminary to show yet, the slope ap- 
pears to be much steeper than the computed gluon slopes 
because hadronization further softens the spectra as dis- 
cussed later. 



C. Transport opacity dependence 

In this Section we provide a qualitative explanation 
for the remarkable invariance of the results on the ac- 
tual angular dependence of the differential cross section 
(2). The reason that this simplification occurs is that for 
the /i and gluon energy range considered in these plots, 
the transport opacity is actually high enough that little 
memory of the initial gluon momentum direction remains 
after multiple scattering. Thus, isotropic scattering and 
/j, = T forward scattering both lead to essentially a ran- 
dom reorientation of all the gluons involved. 

In Fig. 12, aside from a small 10% delta function com- 
ponent due to the gluons in the "corona" surface region 
that escape without rescattering, the bulk of the minijets 
undergo enough rescatterings that their final direction is 
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randomized. In Fig. 13 we also show that the rapid- 
ity shift of gluons in each transverse momentum interval 
considered has the form close to one expected if local 
thermal equilibrium occurred: 

dP f 

= J dm 2 Li dy i dm 2 Lf m± i cosily m±f cosh(y; + Ay) 

x _J J_ e -(m il coshy i )/7 , , e -[m_ L/ cosh(y,+Ay)]/T / 



ATfATf 
Ay cosh Ay — sinh Ay 
sinh 3 Ay 



(14) 



The randomization of momenta is however not suffi- 
cient to ensure the validity of local equilibrium neces- 
sary for the applicability of nondissipative (Euler) hy- 
drodynamics. This is proven by the dependence of the 
transverse momentum spectra and elliptic flow on the fi- 
nite opacity parameter itself. In the hydrodynamic limit 
X = oo and the transport evolution is identical to the hy- 
drodynamical evolution. However, we showed in detail in 
a previous study [17] that the solutions of the transport 
equation still differ very much from ideal hydrodynam- 
ics, for physically extreme x ~ 20 opacities. While no 
covariant 3+1D Navier-Stokes solutions are yet known, 
our transport solutions demonstrate the effects of dissi- 
pation through their dependence on 1/x- 

The invariance of the transport solutions to the angu- 
lar distributions for a fixed \ indicate however that we 
are not extremely far from the local thermal though dis- 
sipativc limit. In particular, our high opacity solutions 
are far from the Eikonal (Knudsen) type dynamics as 
considered in [13]. 

A rough criterion for the validity of the Eikonal ap- 
proximation is that the angle between the initial and 
final parton momenta in the laboratory frame satisfies 
A9 <C 1, say A9 < 0.3. For an energetic parton that 
undergoes N elastic collisions, this angle can be approx- 
imated in an analogous way to random walk as 



(A6) N w yjN{A6 2 ) 

N=l 



(15) 



because the angles transform as A6 w A6 cm \fs/2E. In 
the small angle approximation, we can use Eq. (3) to 
estimate 



(AO 



gt(f) 



Hence, for elastic collisions off typical thermal partons 
(s w GET, n = T), the Eikonal approximation is valid 
if <7 t (s)/<7 < E/(15TN). For N ~ 10, this condition is 
satisfied for E > 20T. Note that the total cross section 
does not depend on the parton energy and therefore the 
number of collisions is approximately independent of en- 
ergy. In Fig. 12 we see that the Eikonal limit is only 
approached slowly as the parton energy increases. 

We conclude from these results that the pattern of "jet 
quenching" , as observed at RHIC via the suppression of 



moderate transverse momentum particles and the satu- 
ration of elliptic flow above some critical pj_ , can be re- 
produced if sufficiently high transport opacities are pos- 
tulated with any angular distribution. 



IV. HADRON SPECTRA 

The results in the previous Section pertain only to par- 
tons. To compare with experimental results, we must 
adopt a model of hadronization. Here we compute the 
hadronic observables from two different hadronization 
schemes. 



A. Hadronization via local parton- hadron duality 

The simplest hadronization scheme is based on the idea 
of local parton- hadron duality. If as in Ref. [21], we as- 
sume that each gluon gets converted to a pion with equal 
probability for the three isospin states, then we may ap- 
proximate the transverse momentum distribution of neg- 
ative charged hadrons roughly as 



A-(p±) 



/tt-(P-L) = »fg{p±)- 



(16) 



With the above prescription, elliptic flow does not 
change during hadronization, i.e., Figs. 4-8 show the 
negative hadron flow as well. Furthermore, the nega- 
tive hadron p± spectra can be obtained from Figs. 9-11 
via simply dividing by 3. Consequently, the scaling (11) 
holds for the negative hadron flow and spectra as well. 

In Fig. 7, the elliptic flow data by STAR are re- 
produced with a transport opacity Xb=a = 47.8. For 
an initial gluon density dN g (0)/dri = 1000, this corre- 
sponds to at ~ 14 mb, i.e., to a total cross section of 
<jq w 45 mb with \x = Tq. If wc took, on the other 
hand, the pQCD gg cross section of 3 mb with \i = Tq, 
this opacity would correspond to an initial gluon density 
of dN g (0)/dn ~ 15000 that is contradicted by the much 
smaller observed dN c h/dn w 600. 

The p± spectra provide a much stronger constraint on 
the initial gluon density as their absolute magnitude is 
proportional to it. At high opacities (x >^ 10), the 
need for high particle subdivisions poses a severe com- 
putational problem, therefore we can reliably compute 
particle spectra for semicentral collisions only. Never- 
theless, the data measured by STAR in central collisions, 
where quenching due to parton energy loss is expected to 
be maximal, provides an important lower bound on the 
particle yields. Figure 11 shows that the elastic transport 
opacities \ <^ 20 considered here are compatible with 
this lower bound. 

On the other hand, the cascade semicentral results are 
very much above the preliminary STAR spectra [30] for 
central collisions. The problem is that the fragmentation 
of quarks and gluons generally soften considerably the 
high-pj_ spectra as we show in the next Section. 
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B. Hadronization via independent fragmentation 

The next simplest hadronization scheme is the frag- 
mentation of gluons as independent jets. Wc consider 
here only the g — > tt channel with the next-to-leading- 
order fragmentation function taken from Ref. [28]. Wc 
took the scale factor s = log(Q 2 )/ log(Qo) to be zero 
because the initial HIJING gluon distribution is already 
"self-quenched" due to initial and final state radiation. 
Also, since we do not consider the contribution of low-pj_ 
soft multiparticle production (beam jet fragmentation), 
we consider from now on only hadrons with p± > 2 GeV. 

Figure 14 shows the final impact-parameter-averaged 
negative hadron flow as a function of the transport opac- 
ity. The flow pattern and the magnitude of the flow 
are much the same as for partons in Fig. 7. Hence, 
we get the same constraint on the initial parameters as 
for hadronization via local parton-hadron duality. In the 
p± < 2 GeV region this simple calculation does not repro- 
duce the data because it does not include contributions 
coming from soft physics. 

The p± spectra of charged hadrons are shown as a func- 
tion of transport opacity and impact parameter in Figs. 
15 and 16. In addition to quenching because of energy 
loss, the final pion spectra are further quenched due to in- 
dependent fragmentation. With this additional quench- 
ing, the parton cascade results approach the preliminary 
STAR data [30], as indicated in Fig. 17, only for rather 
extreme a (n = T) w 100 mb if HIJING dN g /dy = 210 
is assumed, or if a w 25 mb and EKRT dN g /dy = 1000 
is assumed. These elastic cross sections exceed the con- 
ventional few mb pQCD cross sections at this scale by at 
least an order of magnitude. 

V. CONCLUSIONS 

The MPC parton cascade technique was applied to 
solve the covariant Boltzmann transport numerically and 
compute new observables at RHIC. Our focus was on the 
preliminary differential elliptic flow and charged hadron 
moderate p± > 2 GeV/c spectra. We compared results 
using two different hadronization schemes: independent 
fragmentation and local parton-hadron duality. 

Our main result is that if only elastic scattering is 
taken into account in the covariant Boltzmann equa- 
tion, extremely large densities and/or elastic parton 
cross sections, a tot dNjdr\ <~ 80 times the HIJING esti- 
mate, are needed to reproduce the elliptic flow data [16]. 
Hadronization via local parton-hadron duality fails to re- 
produce the rapidly falling high-pj_ spectra. However, in- 
dependent fragmentation of our MPC solutions compare 
well to the charged hadron p± spectra when rather large 
clastic transport opacities are postulated. 

The solutions clearly demonstrate how finite (even ex- 
treme) reaction rates in A + A lead to major deviations 



from ideal hydrodynamic transverse flow effects at trans- 
verse momentum px > 2 GeV. The pattern of quenching 
found with MPC is surprisingly similar to that obtained 
in the two component model of GVW [13]. The main 
difference between the high opacity MPC solutions re- 
ported here and the low opacity results of GVW is that 
the latter include radiative energy loss in an Eikonal for- 
malism joined to a parametrized phenomenological "hy- 
drodynamic" component. 

It is known that radiative energy loss of ultrarelativis- 
tic partons is much larger than elastic energy loss in 
a medium for a fixed cross section. In GVW a simi- 
lar quenching pattern was obtained with more modest 
initial densities dN g /dy ~ 500 and small pQCD elastic 
rates because the induced gluon radiation associated with 
multiple elastic collisions is large enough to compensate 
for the small elastic transport opacity in that case. In 
MPC the same level of quenching is achieved only when 
the elastic opacity is increased artificially by an order of 
magnitude. Therefore, the present study confirms the 
expectation that elastic scattering alone is not enough to 
generate the degree of collectivity observed now at RHIC. 

VI. OUTLOOK 

The results presented here underscore the urgent need 
to develop practical convergent algorithms to incorpo- 
rate inelastic 2 <-» 3 processes. Preliminary work in Ref. 
[18] indicated a rather slow convergence towards Lorentz 
covariance using the particle subdivision technique. Un- 
like the £ -1 / 2 convergence of 2 — > 2 transport solutions, a 
much slower rate of convergence oc is expected with 

the parton subdivision method used to retain Lorentz co- 
variance of 2 <-» 3 processes. 

In addition, a more powerful covariant approximation 
to Boltzmann transport theory may be needed to over- 
come the overwhelming computational difficulty in the 
high opacity regime for central collisions. Wc found that 
even for the case of elastic scattering, particle subdivi- 
sions up to 1000 are required to maintain covariance and 
stabilize the final spectra. 

Finally, we note that all results in this paper pertain 
to slowly varying, smooth initial conditions. In Ref. [29], 
it was suggested that copious minijet production may in- 
duce large (nonstatistical) local fluctuations that could 
evolve in a turbulent manner. A transport study of dif- 
ferential elliptic flow and jet quenching in such inhomoge- 
neous initial conditions would be interesting to compare 
to the presently known hydrodynamic and Boltzmann 
solutions. 
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APPENDIX: COMPARISON OF ELLIPTIC FLOW 
DEFINITIONS 

The main reason that definition (12) gives a very sim- 
ilar elliptic flow to definition (13), contrary to the op- 
posite observation from hydrodynamics, is the difference 
between the hydro and the cascade V2(b) shapes. Hy- 
drodynamical models predict an increasing v 2 out to 
b m 12 — 13 fm, while the cascade v 2 peaks at b m w 8 
fm. The b < b m region where v 2 is small gets a larger 
weight in the STAR definition, while the b > b m region 
where v 2 is as well small gets a smaller weight. The for- 
mer effect tends to reduce elliptic flow, while the latter 
tends to increase it. In the case of hydrodynamics, the 
first effect dominates, while for the cascade, the second 
effect turns out to be larger. 

We illustrate this with a simple analytic calculation. 
The impact parameter dependence of elliptic flow can be 
fitted with the general form 



v 2 {b)=K 



B 



1 



(Al) 



while the particle spectrum is approximately linear in the 
b = 2 — 12 fm region of interest 



dN 
dydp± 



(b) = C 1 



b 

~~B' 



(A2) 



Here the parameters K, B, B', c, d, and C in general 
depend on p±. 

It is easy to show that for the above functions definition 
(12) gives an elliptic flow v% v = 2K T(c+2)T(d+l)/T(c+ 
d + 3), while definition (13) yields v$ TAR = 6K T(c + 
2)T(d + 2)/T(c + d + 4), provided we assume B = B' and 
integrate up to b max = B in both cases. Therefore, 



.STAR 



3d + 3 
c + d + 3 



(A3) 



which is larger than one, if and only if 2d > c. For our 
cascade results, the fits to v 2 (b) give d/c ~ 0.6 — 1.4 
[c(pj_) ~ 1.1 — 3.3, d(p±) <~ 1.0 — 4.5], therefore we have 
v 2 TAR {pi_) > v 2 v (p±). On the other hand, hydrody- 
namics gives an approximately linear v 2 (b), i.e., c = 1 
and d = 0, which results in v 2 (j>±) < v 2 v (p±). 

We found essentially the same when taking an expo- 
nential fit to dN/dp±(b) instead of a linear one. 



Equation (A3) yields v 2 TAR /v 2 v ~ 1.1—1.5 for the cas- 
cade (depending onpi and initial conditions). However, 
in reality B' is smaller than B (typically B' w 11 — 13 
fm, while B w 13— 16 fm), which influences v 2 TAR . Fur- 
thermore, the upper limit of integration b max = 12 fm is 
also smaller than B. This affects primarily v 2 v , through 
the normalization constant l/b ma x- The integrals in both 
definitions are to a large degree insensitive to variations 
of b max because the integrands cut off naturally at large 
b. 

To illustrate these effects, we repeat the previous ana- 
lytic calculation with b max = xB for v 2 v , while B' = yB 
and b m ax = x'B' for v 2 TAR . The integrals yield 



„STAR ( , x _ a j^ y B vx' (c + 2, d + 1) - B yx > (c + 3, d + 1) 
h ( X >y>- bK yV 2 (3-2x') 



v%"(x) = 2K 



B x (c + 2,d+l) 



(A4) 



where B z (a,b) = f^dt t a 1 (l — t) b 1 is the incomplete 
beta function. Thus, 



,STAR 



(l-x',l-y) 



2d 



.STAR 



(1.1) 



l + 1 — T y + 0(x' 2 ) + 0(y 2 



v% v (l - x) 

«r(i) 



= l + 2x + 0(x min { 2 ' d+1 >) 



(A5) 



i.e., the leading correction to v 2 v comes at first order in 
(bmax — B), while the leading correction to v 2 TAR comes 
at first order in (B' — B). 

For parameters and initial conditions considered in this 
study, the corrected formulas (A4) yield v 2 TAR /v 2 v <~ 
0.9 - 1.05 (with x' = 1). The uncertainty of v 2 TAR due 
to the unknown experimental cutoff b ma x (i-e., x') is less 
than 5%. For the hydro (c = 1, d = 0), this analysis gives 
a much smaller ratio v 2 TAR l " %av 
x) k, 0.75. 



M"=3(l-x')(l-y)/4(l- 
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FIG. 1. Strong dependence of the gluon elliptic flow on parton 
subdivision as a function of p ± is shown for Au+Au at y/s = 130A 
GeV with b = 8 fm. Solutions for transport opacity \ = 9-74^ 
(see Table I), and particle subdivisions £ = 1, 5, 50, 225, and 450 
arc shown. 
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FIG. 2. Strong dependence of the final gluon p± spectra on 
parton subdivision is shown for Au+Au at y/s = 130A GeV with 
6 = 8 fm. Solutions for transport opacity \ = 9.74 A ' (see Table I), 
and particle subdivisions I = 1, 5, 50, 225, and 450 are shown as 
in Fig. 1. The spectra are normalized here to dN(Q)/drj = 210. 
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is taken from [6] with the so called sBC initial conditions. It was 
extrapolated beyond p± = 3 GeV using an exponential fit to the 
dN/pj_dpj_ distribution between 2 and 3 GeV. 
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is taken from [6] with the so called sBC initial conditions, and was 
extrapolated beyond p± = 3 GeV using an exponential fit to the 
dN/pj_dpj_ distribution between 2 and 3 GeV. Preliminary STAR 
data [30] suggest that this ratio may reach 10~ 4 in the pj_ ~ 5 — 6 
GeV/c range. 



A) (T = 100 mb, T /n = 1 


B) a = 100 mb, T /n = 


b [ImJ 


(n) 


X 


o [1m] 


(n) 


X 





33.0 


20.2 





35.8 


47.8 


2 


31.7 


19.4 


2 


34.3 


45.8 


4 


28.1 


17.2 


4 


30.2 


40.2 


6 


23.0 


14.1 


6 


24.0 


32.0 


8 


15.9 


9.74 


8 


16.3 


21.8 


10 


8.16 


5.00 


10 


8.23 


11.0 


12 


2.15 


1.32 


12 


2.18 


2.90 


C) do = 40 mb, To/^i = 1 


D) g o = 40 mb, T /n = 


b [fm] 


(n) 


X 


b [fm] 


(n) 


X 





13.4 


8.22 





13.7 


18.3 


2 


12.9 


7.90 


2 


13.2 


17.6 


4 


11.4 


6.98 


4 


11.6 


15.5 


6 


9.26 


5.68 


6 


9.38 


12.5 


8 


6.37 


3.90 


8 


6.44 


8.58 


10 


3.23 


1.98 


10 


3.27 


4.36 


12 


0.86 


0.52 


12 


0.86 


1.14 


E) a = 3 mb, T /(j, = 1 


F) various, b — 8 fm 


b [fm] 


(n) 


X 


(To [fm] 


To In 


(n) 


X 





1.00 


0.62 


60 


1.54 


9.51 


3.68 


2 


0.96 


0.58 


16 





2.55 


3.40 


4 


0.85 


0.52 


100 


1.40 


15.9 


6.86 


6 


0.69 


0.42 


100 


2.21 


15.7 


3.88 


8 


0.47 


0.28 


100 


4.43 


15.5 


1.44 


10 


0.24 


0.148 








12 


0.064 


0.040 









TABLE I. Parameters and transport opacity for each transport 
solution computed via MPC for the present study. 



14 



